<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
  "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">


<html xmlns="http://www.w3.org/1999/xhtml">
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
    
    <title>Chapter 6 Multiple Sequence Alignment objects &mdash; Biopython_en 1.0 documentation</title>
    
    <link rel="stylesheet" href="_static/default.css" type="text/css" />
    <link rel="stylesheet" href="_static/pygments.css" type="text/css" />
    
    <script type="text/javascript">
      var DOCUMENTATION_OPTIONS = {
        URL_ROOT:    './',
        VERSION:     '1.0',
        COLLAPSE_INDEX: false,
        FILE_SUFFIX: '.html',
        HAS_SOURCE:  true
      };
    </script>
    <script type="text/javascript" src="_static/jquery.js"></script>
    <script type="text/javascript" src="_static/underscore.js"></script>
    <script type="text/javascript" src="_static/doctools.js"></script>
    <link rel="top" title="Biopython_en 1.0 documentation" href="index.html" />
    <link rel="next" title="Chapter 7 BLAST" href="chr7.html" />
    <link rel="prev" title="Chapter 5 Sequence Input/Output" href="chr5.html" /> 
  </head>
  <body>
    <div class="related">
      <h3>Navigation</h3>
      <ul>
        <li class="right" style="margin-right: 10px">
          <a href="genindex.html" title="General Index"
             accesskey="I">index</a></li>
        <li class="right" >
          <a href="chr7.html" title="Chapter 7 BLAST"
             accesskey="N">next</a> |</li>
        <li class="right" >
          <a href="chr5.html" title="Chapter 5 Sequence Input/Output"
             accesskey="P">previous</a> |</li>
        <li><a href="index.html">Biopython_en 1.0 documentation</a> &raquo;</li> 
      </ul>
    </div>  

    <div class="document">
      <div class="documentwrapper">
        <div class="bodywrapper">
          <div class="body">
            
  <div class="section" id="chapter-6-multiple-sequence-alignment-objects">
<h1>Chapter 6  Multiple Sequence Alignment objects<a class="headerlink" href="#chapter-6-multiple-sequence-alignment-objects" title="Permalink to this headline">¶</a></h1>
<p>This chapter is about Multiple Sequence Alignments, by which we mean a
collection of multiple sequences which have been aligned together –
usually with the insertion of gap characters, and addition of leading or
trailing gaps – such that all the sequence strings are the same length.
Such an alignment can be regarded as a matrix of letters, where each row
is held as a <tt class="docutils literal"><span class="pre">SeqRecord</span></tt> object internally.</p>
<p>We will introduce the <tt class="docutils literal"><span class="pre">MultipleSeqAlignment</span></tt> object which holds this
kind of data, and the <tt class="docutils literal"><span class="pre">Bio.AlignIO</span></tt> module for reading and writing
them as various file formats (following the design of the <tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt>
module from the previous chapter). Note that both <tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt> and
<tt class="docutils literal"><span class="pre">Bio.AlignIO</span></tt> can read and write sequence alignment files. The
appropriate choice will depend largely on what you want to do with the
data.</p>
<p>The final part of this chapter is about our command line wrappers for
common multiple sequence alignment tools like ClustalW and MUSCLE.</p>
<div class="section" id="parsing-or-reading-sequence-alignments">
<h2>6.1  Parsing or Reading Sequence Alignments<a class="headerlink" href="#parsing-or-reading-sequence-alignments" title="Permalink to this headline">¶</a></h2>
<p>We have two functions for reading in sequence alignments,
<tt class="docutils literal"><span class="pre">Bio.AlignIO.read()</span></tt> and <tt class="docutils literal"><span class="pre">Bio.AlignIO.parse()</span></tt> which following the
convention introduced in <tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt> are for files containing one or
multiple alignments respectively.</p>
<p>Using <tt class="docutils literal"><span class="pre">Bio.AlignIO.parse()</span></tt> will return an <em>iterator</em> which gives
<tt class="docutils literal"><span class="pre">MultipleSeqAlignment</span></tt> objects. Iterators are typically used in a for
loop. Examples of situations where you will have multiple different
alignments include resampled alignments from the PHYLIP tool
<tt class="docutils literal"><span class="pre">seqboot</span></tt>, or multiple pairwise alignments from the EMBOSS tools
<tt class="docutils literal"><span class="pre">water</span></tt> or <tt class="docutils literal"><span class="pre">needle</span></tt>, or Bill Pearson’s FASTA tools.</p>
<p>However, in many situations you will be dealing with files which contain
only a single alignment. In this case, you should use the
<tt class="docutils literal"><span class="pre">Bio.AlignIO.read()</span></tt> function which returns a single
<tt class="docutils literal"><span class="pre">MultipleSeqAlignment</span></tt> object.</p>
<p>Both functions expect two mandatory arguments:</p>
<ol class="arabic simple">
<li>The first argument is a <em>handle</em> to read the data from, typically an
open file (see Section <a class="reference external" href="#sec:appendix-handles">22.1</a>), or a
filename.</li>
<li>The second argument is a lower case string specifying the alignment
format. As in <tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt> we don’t try and guess the file format
for you! See
<tt class="docutils literal"><span class="pre">`http://biopython.org/wiki/AlignIO</span></tt> &lt;<a class="reference external" href="http://biopython.org/wiki/AlignIO">http://biopython.org/wiki/AlignIO</a>&gt;`__
for a full listing of supported formats.</li>
</ol>
<p>There is also an optional <tt class="docutils literal"><span class="pre">seq_count</span></tt> argument which is discussed in
Section <a class="reference external" href="#sec:AlignIO-count-argument">6.1.3</a> below for dealing with
ambiguous file formats which may contain more than one alignment.</p>
<p>A further optional <tt class="docutils literal"><span class="pre">alphabet</span></tt> argument allowing you to specify the
expected alphabet. This can be useful as many alignment file formats do
not explicitly label the sequences as RNA, DNA or protein – which means
<tt class="docutils literal"><span class="pre">Bio.AlignIO</span></tt> will default to using a generic alphabet.</p>
<div class="section" id="single-alignments">
<h3>6.1.1  Single Alignments<a class="headerlink" href="#single-alignments" title="Permalink to this headline">¶</a></h3>
<p>As an example, consider the following annotation rich protein alignment
in the PFAM or Stockholm file format:</p>
<p>This is the seed alignment for the Phage_Coat_Gp8 (PF05371) PFAM
entry, downloaded from a now out of date release of PFAM from
<tt class="docutils literal"><span class="pre">`http://pfam.sanger.ac.uk/</span></tt> &lt;<a class="reference external" href="http://pfam.sanger.ac.uk/">http://pfam.sanger.ac.uk/</a>&gt;`__. We can
load this file as follows (assuming it has been saved to disk as
“PF05371_seed.sth” in the current working directory):</p>
<p>This code will print out a summary of the alignment:</p>
<p>You’ll notice in the above output the sequences have been truncated. We
could instead write our own code to format this as we please by
iterating over the rows as <tt class="docutils literal"><span class="pre">SeqRecord</span></tt> objects:</p>
<p>You could also use the alignment object’s <tt class="docutils literal"><span class="pre">format</span></tt> method to show it
in a particular file format – see
Section <a class="reference external" href="#sec:alignment-format-method">6.2.2</a> for details.</p>
<p>Did you notice in the raw file above that several of the sequences
include database cross-references to the PDB and the associated known
secondary structure? Try this:</p>
<p>To have a look at all the sequence annotation, try this:</p>
<p>Sanger provide a nice web interface at
<tt class="docutils literal"><span class="pre">`http://pfam.sanger.ac.uk/family?acc=PF05371</span></tt> &lt;<a class="reference external" href="http://pfam.sanger.ac.uk/family?acc=PF05371">http://pfam.sanger.ac.uk/family?acc=PF05371</a>&gt;`__
which will actually let you download this alignment in several other
formats. This is what the file looks like in the FASTA file format:</p>
<p>Note the website should have an option about showing gaps as periods
(dots) or dashes, we’ve shown dashes above. Assuming you download and
save this as file “PF05371_seed.faa” then you can load it with almost
exactly the same code:</p>
<p>All that has changed in this code is the filename and the format string.
You’ll get the same output as before, the sequences and record
identifiers are the same. However, as you should expect, if you check
each <tt class="docutils literal"><span class="pre">SeqRecord</span></tt> there is no annotation nor database cross-references
because these are not included in the FASTA file format.</p>
<p>Note that rather than using the Sanger website, you could have used
<tt class="docutils literal"><span class="pre">Bio.AlignIO</span></tt> to convert the original Stockholm format file into a
FASTA file yourself (see below).</p>
<p>With any supported file format, you can load an alignment in exactly the
same way just by changing the format string. For example, use “phylip”
for PHYLIP files, “nexus” for NEXUS files or “emboss” for the alignments
output by the EMBOSS tools. There is a full listing on the wiki page
(<tt class="docutils literal"><span class="pre">`http://biopython.org/wiki/AlignIO</span></tt> &lt;<a class="reference external" href="http://biopython.org/wiki/AlignIO">http://biopython.org/wiki/AlignIO</a>&gt;`__)
and in the built in documentation (also
<a class="reference external" href="http://biopython.org/DIST/docs/api/Bio.AlignIO-module.html">online</a>):</p>
</div>
<div class="section" id="multiple-alignments">
<h3>6.1.2  Multiple Alignments<a class="headerlink" href="#multiple-alignments" title="Permalink to this headline">¶</a></h3>
<p>The previous section focused on reading files containing a single
alignment. In general however, files can contain more than one
alignment, and to read these files we must use the
<tt class="docutils literal"><span class="pre">Bio.AlignIO.parse()</span></tt> function.</p>
<p>Suppose you have a small alignment in PHYLIP format:</p>
<p>If you wanted to bootstrap a phylogenetic tree using the PHYLIP tools,
one of the steps would be to create a set of many resampled alignments
using the tool <tt class="docutils literal"><span class="pre">bootseq</span></tt>. This would give output something like this,
which has been abbreviated for conciseness:</p>
<p>If you wanted to read this in using <tt class="docutils literal"><span class="pre">Bio.AlignIO</span></tt> you could use:</p>
<p>This would give the following output, again abbreviated for display:</p>
<p>As with the function <tt class="docutils literal"><span class="pre">Bio.SeqIO.parse()</span></tt>, using
<tt class="docutils literal"><span class="pre">Bio.AlignIO.parse()</span></tt> returns an iterator. If you want to keep all the
alignments in memory at once, which will allow you to access them in any
order, then turn the iterator into a list:</p>
</div>
<div class="section" id="ambiguous-alignments">
<h3>6.1.3  Ambiguous Alignments<a class="headerlink" href="#ambiguous-alignments" title="Permalink to this headline">¶</a></h3>
<p>Many alignment file formats can explicitly store more than one
alignment, and the division between each alignment is clear. However,
when a general sequence file format has been used there is no such block
structure. The most common such situation is when alignments have been
saved in the FASTA file format. For example consider the following:</p>
<p>This could be a single alignment containing six sequences (with repeated
identifiers). Or, judging from the identifiers, this is probably two
different alignments each with three sequences, which happen to all have
the same length.</p>
<p>What about this next example?</p>
<p>Again, this could be a single alignment with six sequences. However this
time based on the identifiers we might guess this is three pairwise
alignments which by chance have all got the same lengths.</p>
<p>This final example is similar:</p>
<p>In this third example, because of the differing lengths, this cannot be
treated as a single alignment containing all six records. However, it
could be three pairwise alignments.</p>
<p>Clearly trying to store more than one alignment in a FASTA file is not
ideal. However, if you are forced to deal with these as input files
<tt class="docutils literal"><span class="pre">Bio.AlignIO</span></tt> can cope with the most common situation where all the
alignments have the same number of records. One example of this is a
collection of pairwise alignments, which can be produced by the EMBOSS
tools <tt class="docutils literal"><span class="pre">needle</span></tt> and <tt class="docutils literal"><span class="pre">water</span></tt> – although in this situation,
<tt class="docutils literal"><span class="pre">Bio.AlignIO</span></tt> should be able to understand their native output using
“emboss” as the format string.</p>
<p>To interpret these FASTA examples as several separate alignments, we can
use <tt class="docutils literal"><span class="pre">Bio.AlignIO.parse()</span></tt> with the optional <tt class="docutils literal"><span class="pre">seq_count</span></tt> argument
which specifies how many sequences are expected in each alignment (in
these examples, 3, 2 and 2 respectively). For example, using the third
example as the input data:</p>
<p>giving:</p>
<p>Using <tt class="docutils literal"><span class="pre">Bio.AlignIO.read()</span></tt> or <tt class="docutils literal"><span class="pre">Bio.AlignIO.parse()</span></tt> without the
<tt class="docutils literal"><span class="pre">seq_count</span></tt> argument would give a single alignment containing all six
records for the first two examples. For the third example, an exception
would be raised because the lengths differ preventing them being turned
into a single alignment.</p>
<p>If the file format itself has a block structure allowing <tt class="docutils literal"><span class="pre">Bio.AlignIO</span></tt>
to determine the number of sequences in each alignment directly, then
the <tt class="docutils literal"><span class="pre">seq_count</span></tt> argument is not needed. If it is supplied, and doesn’t
agree with the file contents, an error is raised.</p>
<p>Note that this optional <tt class="docutils literal"><span class="pre">seq_count</span></tt> argument assumes each alignment in
the file has the same number of sequences. Hypothetically you may come
across stranger situations, for example a FASTA file containing several
alignments each with a different number of sequences – although I would
love to hear of a real world example of this. Assuming you cannot get
the data in a nicer file format, there is no straight forward way to
deal with this using <tt class="docutils literal"><span class="pre">Bio.AlignIO</span></tt>. In this case, you could consider
reading in the sequences themselves using <tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt> and batching
them together to create the alignments as appropriate.</p>
</div>
</div>
<div class="section" id="writing-alignments">
<h2>6.2  Writing Alignments<a class="headerlink" href="#writing-alignments" title="Permalink to this headline">¶</a></h2>
<p>We’ve talked about using <tt class="docutils literal"><span class="pre">Bio.AlignIO.read()</span></tt> and
<tt class="docutils literal"><span class="pre">Bio.AlignIO.parse()</span></tt> for alignment input (reading files), and now
we’ll look at <tt class="docutils literal"><span class="pre">Bio.AlignIO.write()</span></tt> which is for alignment output
(writing files). This is a function taking three arguments: some
<tt class="docutils literal"><span class="pre">MultipleSeqAlignment</span></tt> objects (or for backwards compatibility the
obsolete <tt class="docutils literal"><span class="pre">Alignment</span></tt> objects), a handle or filename to write to, and a
sequence format.</p>
<p>Here is an example, where we start by creating a few
<tt class="docutils literal"><span class="pre">MultipleSeqAlignment</span></tt> objects the hard way (by hand, rather than by
loading them from a file). Note we create some <tt class="docutils literal"><span class="pre">SeqRecord</span></tt> objects to
construct the alignment from.</p>
<p>Now we have a list of <tt class="docutils literal"><span class="pre">Alignment</span></tt> objects, we’ll write them to a
PHYLIP format file:</p>
<p>And if you open this file in your favourite text editor it should look
like this:</p>
<p>Its more common to want to load an existing alignment, and save that,
perhaps after some simple manipulation like removing certain rows or
columns.</p>
<p>Suppose you wanted to know how many alignments the
<tt class="docutils literal"><span class="pre">Bio.AlignIO.write()</span></tt> function wrote to the handle? If your alignments
were in a list like the example above, you could just use
<tt class="docutils literal"><span class="pre">len(my_alignments)</span></tt>, however you can’t do that when your records come
from a generator/iterator. Therefore the <tt class="docutils literal"><span class="pre">Bio.AlignIO.write()</span></tt>
function returns the number of alignments written to the file.</p>
<p><em>Note</em> - If you tell the <tt class="docutils literal"><span class="pre">Bio.AlignIO.write()</span></tt> function to write to a
file that already exists, the old file will be overwritten without any
warning.</p>
<div class="section" id="converting-between-sequence-alignment-file-formats">
<h3>6.2.1  Converting between sequence alignment file formats<a class="headerlink" href="#converting-between-sequence-alignment-file-formats" title="Permalink to this headline">¶</a></h3>
<p>Converting between sequence alignment file formats with <tt class="docutils literal"><span class="pre">Bio.AlignIO</span></tt>
works in the same way as converting between sequence file formats with
<tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt> (Section <a class="reference external" href="#sec:SeqIO-conversion">5.5.2</a>). We load
generally the alignment(s) using <tt class="docutils literal"><span class="pre">Bio.AlignIO.parse()</span></tt> and then save
them using the <tt class="docutils literal"><span class="pre">Bio.AlignIO.write()</span></tt> – or just use the
<tt class="docutils literal"><span class="pre">Bio.AlignIO.convert()</span></tt> helper function.</p>
<p>For this example, we’ll load the PFAM/Stockholm format file used earlier
and save it as a Clustal W format file:</p>
<p>Or, using <tt class="docutils literal"><span class="pre">Bio.AlignIO.parse()</span></tt> and <tt class="docutils literal"><span class="pre">Bio.AlignIO.write()</span></tt>:</p>
<p>The <tt class="docutils literal"><span class="pre">Bio.AlignIO.write()</span></tt> function expects to be given multiple
alignment objects. In the example above we gave it the alignment
iterator returned by <tt class="docutils literal"><span class="pre">Bio.AlignIO.parse()</span></tt>.</p>
<p>In this case, we know there is only one alignment in the file so we
could have used <tt class="docutils literal"><span class="pre">Bio.AlignIO.read()</span></tt> instead, but notice we have to
pass this alignment to <tt class="docutils literal"><span class="pre">Bio.AlignIO.write()</span></tt> as a single element list:</p>
<p>Either way, you should end up with the same new Clustal W format file
“PF05371_seed.aln” with the following content:</p>
<p>Alternatively, you could make a PHYLIP format file which we’ll name
“PF05371_seed.phy”:</p>
<p>This time the output looks like this:</p>
<p>One of the big handicaps of the PHYLIP alignment file format is that the
sequence identifiers are strictly truncated at ten characters. In this
example, as you can see the resulting names are still unique - but they
are not very readable. In this particular case, there is no clear way to
compress the identifiers, but for the sake of argument you may want to
assign your own names or numbering system. This following bit of code
manipulates the record identifiers before saving the output:</p>
<p>This code used a Python dictionary to record a simple mapping from the
new sequence system to the original identifier:</p>
<p>Here is the new PHYLIP format output:</p>
<p>In general, because of the identifier limitation, working with PHYLIP
file formats shouldn’t be your first choice. Using the PFAM/Stockholm
format on the other hand allows you to record a lot of additional
annotation too.</p>
</div>
<div class="section" id="getting-your-alignment-objects-as-formatted-strings">
<h3>6.2.2  Getting your alignment objects as formatted strings<a class="headerlink" href="#getting-your-alignment-objects-as-formatted-strings" title="Permalink to this headline">¶</a></h3>
<p>The <tt class="docutils literal"><span class="pre">Bio.AlignIO</span></tt> interface is based on handles, which means if you
want to get your alignment(s) into a string in a particular file format
you need to do a little bit more work (see below). However, you will
probably prefer to take advantage of the alignment object’s <tt class="docutils literal"><span class="pre">format()</span></tt>
method. This takes a single mandatory argument, a lower case string
which is supported by <tt class="docutils literal"><span class="pre">Bio.AlignIO</span></tt> as an output format. For example:</p>
<p>As described in Section <a class="reference external" href="#sec:SeqRecord-format">4.5</a>), the
<tt class="docutils literal"><span class="pre">SeqRecord</span></tt> object has a similar method using output formats supported
by <tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt>.</p>
<p>Internally the <tt class="docutils literal"><span class="pre">format()</span></tt> method is using the <tt class="docutils literal"><span class="pre">StringIO</span></tt> string
based handle and calling <tt class="docutils literal"><span class="pre">Bio.AlignIO.write()</span></tt>. You can do this in
your own code if for example you are using an older version of
Biopython:</p>
</div>
</div>
<div class="section" id="manipulating-alignments">
<h2>6.3  Manipulating Alignments<a class="headerlink" href="#manipulating-alignments" title="Permalink to this headline">¶</a></h2>
<p>Now that we’ve covered loading and saving alignments, we’ll look at what
else you can do with them.</p>
<div class="section" id="slicing-alignments">
<h3>6.3.1  Slicing alignments<a class="headerlink" href="#slicing-alignments" title="Permalink to this headline">¶</a></h3>
<p>First of all, in some senses the alignment objects act like a Python
<tt class="docutils literal"><span class="pre">list</span></tt> of <tt class="docutils literal"><span class="pre">SeqRecord</span></tt> objects (the rows). With this model in mind
hopefully the actions of <tt class="docutils literal"><span class="pre">len()</span></tt> (the number of rows) and iteration
(each row as a <tt class="docutils literal"><span class="pre">SeqRecord</span></tt>) make sense:</p>
<p>You can also use the list-like <tt class="docutils literal"><span class="pre">append</span></tt> and <tt class="docutils literal"><span class="pre">extend</span></tt> methods to add
more rows to the alignment (as <tt class="docutils literal"><span class="pre">SeqRecord</span></tt> objects). Keeping the list
metaphor in mind, simple slicing of the alignment should also make sense
- it selects some of the rows giving back another alignment object:</p>
<p>What if you wanted to select by column? Those of you who have used the
NumPy matrix or array objects won’t be surprised at this - you use a
double index.</p>
<p>Using two integer indices pulls out a single letter, short hand for
this:</p>
<p>You can pull out a single column as a string like this:</p>
<p>You can also select a range of columns. For example, to pick out those
same three rows we extracted earlier, but take just their first six
columns:</p>
<p>Leaving the first index as <tt class="docutils literal"><span class="pre">:</span></tt> means take all the rows:</p>
<p>This brings us to a neat way to remove a section. Notice columns 7, 8
and 9 which are gaps in three of the seven sequences:</p>
<p>Again, you can slice to get everything after the ninth column:</p>
<p>Now, the interesting thing is that addition of alignment objects works
by column. This lets you do this as a way to remove a block of columns:</p>
<p>Another common use of alignment addition would be to combine alignments
for several different genes into a meta-alignment. Watch out though -
the identifiers need to match up (see
Section <a class="reference external" href="#sec:SeqRecord-addition">4.7</a> for how adding <tt class="docutils literal"><span class="pre">SeqRecord</span></tt>
objects works). You may find it helpful to first sort the alignment rows
alphabetically by id:</p>
<p>Note that you can only add two alignments together if they have the same
number of rows.</p>
</div>
<div class="section" id="alignments-as-arrays">
<h3>6.3.2  Alignments as arrays<a class="headerlink" href="#alignments-as-arrays" title="Permalink to this headline">¶</a></h3>
<p>Depending on what you are doing, it can be more useful to turn the
alignment object into an array of letters – and you can do this with
NumPy:</p>
<p>If you will be working heavily with the columns, you can tell NumPy to
store the array by column (as in Fortran) rather then its default of by
row (as in C):</p>
<p>Note that this leaves the original Biopython alignment object and the
NumPy array in memory as separate objects - editing one will not update
the other!</p>
</div>
</div>
<div class="section" id="alignment-tools">
<h2>6.4  Alignment Tools<a class="headerlink" href="#alignment-tools" title="Permalink to this headline">¶</a></h2>
<p>There are <em>lots</em> of algorithms out there for aligning sequences, both
pairwise alignments and multiple sequence alignments. These calculations
are relatively slow, and you generally wouldn’t want to write such an
algorithm in Python. Instead, you can use Biopython to invoke a command
line tool on your behalf. Normally you would:</p>
<ol class="arabic simple">
<li>Prepare an input file of your unaligned sequences, typically this
will be a FASTA file which you might create using <tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt> (see
Chapter <a class="reference external" href="#chapter:Bio.SeqIO">5</a>).</li>
<li>Call the command line tool to process this input file, typically via
one of Biopython’s command line wrappers (which we’ll discuss here).</li>
<li>Read the output from the tool, i.e. your aligned sequences, typically
using <tt class="docutils literal"><span class="pre">Bio.AlignIO</span></tt> (see earlier in this chapter).</li>
</ol>
<p>All the command line wrappers we’re going to talk about in this chapter
follow the same style. You create a command line object specifying the
options (e.g. the input filename and the output filename), then invoke
this command line via a Python operating system call (e.g. using the
<tt class="docutils literal"><span class="pre">subprocess</span></tt> module).</p>
<p>Most of these wrappers are defined in the <tt class="docutils literal"><span class="pre">Bio.Align.Applications</span></tt>
module:</p>
<p>(Ignore the entries starting with an underscore – these have special
meaning in Python.) The module <tt class="docutils literal"><span class="pre">Bio.Emboss.Applications</span></tt> has wrappers
for some of the <a class="reference external" href="http://emboss.sourceforge.net/">EMBOSS suite</a>,
including <tt class="docutils literal"><span class="pre">needle</span></tt> and <tt class="docutils literal"><span class="pre">water</span></tt>, which are described below in
Section <a class="reference external" href="#seq:emboss-needle-water">6.4.5</a>, and wrappers for the
EMBOSS packaged versions of the PHYLIP tools (which EMBOSS refer to as
one of their EMBASSY packages - third party tools with an EMBOSS style
interface). We won’t explore all these alignment tools here in the
section, just a sample, but the same principles apply.</p>
<div class="section" id="clustalw">
<h3>6.4.1  ClustalW<a class="headerlink" href="#clustalw" title="Permalink to this headline">¶</a></h3>
<p>ClustalW is a popular command line tool for multiple sequence alignment
(there is also a graphical interface called ClustalX). Biopython’s
<tt class="docutils literal"><span class="pre">Bio.Align.Applications</span></tt> module has a wrapper for this alignment tool
(and several others).</p>
<p>Before trying to use ClustalW from within Python, you should first try
running the ClustalW tool yourself by hand at the command line, to
familiarise yourself the other options. You’ll find the Biopython
wrapper is very faithful to the actual command line API:</p>
<p>For the most basic usage, all you need is to have a FASTA input file,
such as
<a class="reference external" href="http://biopython.org/DIST/docs/tutorial/examples/opuntia.fasta">opuntia.fasta</a>
(available online or in the Doc/examples subdirectory of the Biopython
source code). This is a small FASTA file containing seven prickly-pear
DNA sequences (from the cactus family <em>Opuntia</em>).</p>
<p>By default ClustalW will generate an alignment and guide tree file with
names based on the input FASTA file, in this case <tt class="docutils literal"><span class="pre">opuntia.aln</span></tt> and
<tt class="docutils literal"><span class="pre">opuntia.dnd</span></tt>, but you can override this or make it explicit:</p>
<p>Notice here we have given the executable name as <tt class="docutils literal"><span class="pre">clustalw2</span></tt>,
indicating we have version two installed, which has a different filename
to version one (<tt class="docutils literal"><span class="pre">clustalw</span></tt>, the default). Fortunately both versions
support the same set of arguments at the command line (and indeed,
should be functionally identical).</p>
<p>You may find that even though you have ClustalW installed, the above
command doesn’t work – you may get a message about “command not found”
(especially on Windows). This indicated that the ClustalW executable is
not on your PATH (an environment variable, a list of directories to be
searched). You can either update your PATH setting to include the
location of your copy of ClustalW tools (how you do this will depend on
your OS), or simply type in the full path of the tool. For example:</p>
<p>Remember, in Python strings <tt class="docutils literal"><span class="pre">\n</span></tt> and <tt class="docutils literal"><span class="pre">\t</span></tt> are by default interpreted
as a new line and a tab – which is why we’re put a letter “r” at the
start for a raw string that isn’t translated in this way. This is
generally good practice when specifying a Windows style file name.</p>
<p>Internally this uses the <tt class="docutils literal"><span class="pre">subprocess</span></tt> module which is now the
recommended way to run another program in Python. This replaces older
options like the <tt class="docutils literal"><span class="pre">os.system()</span></tt> and the <tt class="docutils literal"><span class="pre">os.popen*</span></tt> functions.</p>
<p>Now, at this point it helps to know about how command line tools “work”.
When you run a tool at the command line, it will often print text output
directly to screen. This text can be captured or redirected, via two
“pipes”, called standard output (the normal results) and standard error
(for error messages and debug messages). There is also standard input,
which is any text fed into the tool. These names get shortened to stdin,
stdout and stderr. When the tool finishes, it has a return code (an
integer), which by convention is zero for success.</p>
<p>When you run the command line tool like this via the Biopython wrapper,
it will wait for it to finish, and check the return code. If this is non
zero (indicating an error), an exception is raised. The wrapper then
returns two strings, stdout and stderr.</p>
<p>In the case of ClustalW, when run at the command line all the important
output is written directly to the output files. Everything normally
printed to screen while you wait (via stdout or stderr) is boring and
can be ignored (assuming it worked).</p>
<p>What we care about are the two output files, the alignment and the guide
tree. We didn’t tell ClustalW what filenames to use, but it defaults to
picking names based on the input file. In this case the output should be
in the file <tt class="docutils literal"><span class="pre">opuntia.aln</span></tt>. You should be able to work out how to read
in the alignment using <tt class="docutils literal"><span class="pre">Bio.AlignIO</span></tt> by now:</p>
<p>In case you are interested (and this is an aside from the main thrust of
this chapter), the <tt class="docutils literal"><span class="pre">opuntia.dnd</span></tt> file ClustalW creates is just a
standard Newick tree file, and <tt class="docutils literal"><span class="pre">Bio.Phylo</span></tt> can parse these:</p>
<p>Chapter <a class="reference external" href="#sec:Phylo">13</a> covers Biopython’s support for phylogenetic
trees in more depth.</p>
</div>
<div class="section" id="muscle">
<h3>6.4.2  MUSCLE<a class="headerlink" href="#muscle" title="Permalink to this headline">¶</a></h3>
<p>MUSCLE is a more recent multiple sequence alignment tool than ClustalW,
and Biopython also has a wrapper for it under the
<tt class="docutils literal"><span class="pre">Bio.Align.Applications</span></tt> module. As before, we recommend you try using
MUSCLE from the command line before trying it from within Python, as the
Biopython wrapper is very faithful to the actual command line API:</p>
<p>For the most basic usage, all you need is to have a FASTA input file,
such as
<a class="reference external" href="http://biopython.org/DIST/docs/tutorial/examples/opuntia.fasta">opuntia.fasta</a>
(available online or in the Doc/examples subdirectory of the Biopython
source code). You can then tell MUSCLE to read in this FASTA file, and
write the alignment to an output file:</p>
<p>Note that MUSCLE uses “-in” and “-out” but in Biopython we have to use
“input” and “out” as the keyword arguments or property names. This is
because “in” is a reserved word in Python.</p>
<p>By default MUSCLE will output the alignment as a FASTA file (using
gapped sequences). The <tt class="docutils literal"><span class="pre">Bio.AlignIO</span></tt> module should be able to read
this alignment using <tt class="docutils literal"><span class="pre">format=&quot;fasta&quot;</span></tt>. You can also ask for
ClustalW-like output:</p>
<p>Or, strict ClustalW output where the original ClustalW header line is
used for maximum compatibility:</p>
<p>The <tt class="docutils literal"><span class="pre">Bio.AlignIO</span></tt> module should be able to read these alignments using
<tt class="docutils literal"><span class="pre">format=&quot;clustal&quot;</span></tt>.</p>
<p>MUSCLE can also output in GCG MSF format (using the <tt class="docutils literal"><span class="pre">msf</span></tt> argument),
but Biopython can’t currently parse that, or using HTML which would give
a human readable web page (not suitable for parsing).</p>
<p>You can also set the other optional parameters, for example the maximum
number of iterations. See the built in help for details.</p>
<p>You would then run MUSCLE command line string as described above for
ClustalW, and parse the output using <tt class="docutils literal"><span class="pre">Bio.AlignIO</span></tt> to get an alignment
object.</p>
</div>
<div class="section" id="muscle-using-stdout">
<h3>6.4.3  MUSCLE using stdout<a class="headerlink" href="#muscle-using-stdout" title="Permalink to this headline">¶</a></h3>
<p>Using a MUSCLE command line as in the examples above will write the
alignment to a file. This means there will be no important information
written to the standard out (stdout) or standard error (stderr) handles.
However, by default MUSCLE will write the alignment to standard output
(stdout). We can take advantage of this to avoid having a temporary
output file! For example:</p>
<p>If we run this via the wrapper, we get back the output as a string. In
order to parse this we can use <tt class="docutils literal"><span class="pre">StringIO</span></tt> to turn it into a handle.
Remember that MUSCLE defaults to using FASTA as the output format:</p>
<p>The above approach is fairly simple, but if you are dealing with very
large output text the fact that all of stdout and stderr is loaded into
memory as a string can be a potential drawback. Using the <tt class="docutils literal"><span class="pre">subprocess</span></tt>
module we can work directly with handles instead:</p>
</div>
<div class="section" id="muscle-using-stdin-and-stdout">
<h3>6.4.4  MUSCLE using stdin and stdout<a class="headerlink" href="#muscle-using-stdin-and-stdout" title="Permalink to this headline">¶</a></h3>
<p>We don’t actually <em>need</em> to have our FASTA input sequences prepared in a
file, because by default MUSCLE will read in the input sequence from
standard input! Note this is a bit more advanced and fiddly, so don’t
bother with this technique unless you need to.</p>
<p>First, we’ll need some unaligned sequences in memory as <tt class="docutils literal"><span class="pre">SeqRecord</span></tt>
objects. For this demonstration I’m going to use a filtered version of
the original FASTA file (using a generator expression), taking just six
of the seven sequences:</p>
<p>Then we create the MUSCLE command line, leaving the input and output to
their defaults (stdin and stdout). I’m also going to ask for strict
ClustalW format as for the output.</p>
<p>Now for the fiddly bits using the <tt class="docutils literal"><span class="pre">subprocess</span></tt> module, stdin and
stdout:</p>
<p>That should start MUSCLE, but it will be sitting waiting for its FASTA
input sequences, which we must supply via its stdin handle:</p>
<p>After writing the six sequences to the handle, MUSCLE will still be
waiting to see if that is all the FASTA sequences or not – so we must
signal that this is all the input data by closing the handle. At that
point MUSCLE should start to run, and we can ask for the output:</p>
<p>Wow! There we are with a new alignment of just the six records, without
having created a temporary FASTA input file, or a temporary alignment
output file. However, a word of caution: Dealing with errors with this
style of calling external programs is much more complicated. It also
becomes far harder to diagnose problems, because you can’t try running
MUSCLE manually outside of Biopython (because you don’t have the input
file to supply). There can also be subtle cross platform issues (e.g.
Windows versus Linux), and how you run your script can have an impact
(e.g. at the command line, from IDLE or an IDE, or as a GUI script).
These are all generic Python issues though, and not specific to
Biopython.</p>
<p>If you find working directly with <tt class="docutils literal"><span class="pre">subprocess</span></tt> like this scary, there
is an alternative. If you execute the tool with <tt class="docutils literal"><span class="pre">muscle_cline()</span></tt> you
can supply any standard input as a big string,
<tt class="docutils literal"><span class="pre">muscle_cline(stdin=...)</span></tt>. So, provided your data isn’t very big, you
can prepare the FASTA input in memory as a string using <tt class="docutils literal"><span class="pre">StringIO</span></tt>
(see Section <a class="reference external" href="#sec:appendix-handles">22.1</a>):</p>
<p>You can then run the tool and parse the alignment as follows:</p>
<p>You might find this easier, but it does require more memory (RAM) for
the strings used for the input FASTA and output Clustal formatted data.</p>
</div>
<div class="section" id="emboss-needle-and-water">
<h3>6.4.5  EMBOSS needle and water<a class="headerlink" href="#emboss-needle-and-water" title="Permalink to this headline">¶</a></h3>
<p>The <a class="reference external" href="http://emboss.sourceforge.net/">EMBOSS</a> suite includes the
<tt class="docutils literal"><span class="pre">water</span></tt> and <tt class="docutils literal"><span class="pre">needle</span></tt> tools for Smith-Waterman algorithm local
alignment, and Needleman-Wunsch global alignment. The tools share the
same style interface, so switching between the two is trivial – we’ll
just use <tt class="docutils literal"><span class="pre">needle</span></tt> here.</p>
<p>Suppose you want to do a global pairwise alignment between two
sequences, prepared in FASTA format as follows:</p>
<p>in a file <tt class="docutils literal"><span class="pre">alpha.fasta</span></tt>, and secondly in a file <tt class="docutils literal"><span class="pre">beta.fasta</span></tt>:</p>
<p>Let’s start by creating a complete <tt class="docutils literal"><span class="pre">needle</span></tt> command line object in one
go:</p>
<p>Why not try running this by hand at the command prompt? You should see
it does a pairwise comparison and records the output in the file
<tt class="docutils literal"><span class="pre">needle.txt</span></tt> (in the default EMBOSS alignment file format).</p>
<p>Even if you have EMBOSS installed, running this command may not work –
you might get a message about “command not found” (especially on
Windows). This probably means that the EMBOSS tools are not on your PATH
environment variable. You can either update your PATH setting, or simply
tell Biopython the full path to the tool, for example:</p>
<p>Remember in Python that for a default string <tt class="docutils literal"><span class="pre">\n</span></tt> or <tt class="docutils literal"><span class="pre">\t</span></tt> means a
new line or a tab – which is why we’re put a letter “r” at the start for
a raw string.</p>
<p>At this point it might help to try running the EMBOSS tools yourself by
hand at the command line, to familiarise yourself the other options and
compare them to the Biopython help text:</p>
<p>Note that you can also specify (or change or look at) the settings like
this:</p>
<p>Next we want to use Python to run this command for us. As explained
above, for full control, we recommend you use the built in Python
<tt class="docutils literal"><span class="pre">subprocess</span></tt> module, but for simple usage the wrapper object usually
suffices:</p>
<p>Next we can load the output file with <tt class="docutils literal"><span class="pre">Bio.AlignIO</span></tt> as discussed
earlier in this chapter, as the <tt class="docutils literal"><span class="pre">emboss</span></tt> format:</p>
<p>In this example, we told EMBOSS to write the output to a file, but you
<em>can</em> tell it to write the output to stdout instead (useful if you don’t
want a temporary output file to get rid of – use <tt class="docutils literal"><span class="pre">stdout=True</span></tt> rather
than the <tt class="docutils literal"><span class="pre">outfile</span></tt> argument), and also to read <em>one</em> of the one of the
inputs from stdin (e.g. <tt class="docutils literal"><span class="pre">asequence=&quot;stdin&quot;</span></tt>, much like in the MUSCLE
example in the section above).</p>
<p>This has only scratched the surface of what you can do with <tt class="docutils literal"><span class="pre">needle</span></tt>
and <tt class="docutils literal"><span class="pre">water</span></tt>. One useful trick is that the second file can contain
multiple sequences (say five), and then EMBOSS will do five pairwise
alignments.</p>
<p>Note - Biopython includes its own pairwise alignment code in the
<tt class="docutils literal"><span class="pre">Bio.pairwise2</span></tt> module (written in C for speed, but with a pure Python
fallback available too). This doesn’t work with alignment objects, so we
have not covered it within this chapter. See the module’s docstring
(built in help) for details.</p>
</div>
</div>
</div>


          </div>
        </div>
      </div>
      <div class="sphinxsidebar">
        <div class="sphinxsidebarwrapper">
  <h3><a href="index.html">Table Of Contents</a></h3>
  <ul>
<li><a class="reference internal" href="#">Chapter 6  Multiple Sequence Alignment objects</a><ul>
<li><a class="reference internal" href="#parsing-or-reading-sequence-alignments">6.1  Parsing or Reading Sequence Alignments</a><ul>
<li><a class="reference internal" href="#single-alignments">6.1.1  Single Alignments</a></li>
<li><a class="reference internal" href="#multiple-alignments">6.1.2  Multiple Alignments</a></li>
<li><a class="reference internal" href="#ambiguous-alignments">6.1.3  Ambiguous Alignments</a></li>
</ul>
</li>
<li><a class="reference internal" href="#writing-alignments">6.2  Writing Alignments</a><ul>
<li><a class="reference internal" href="#converting-between-sequence-alignment-file-formats">6.2.1  Converting between sequence alignment file formats</a></li>
<li><a class="reference internal" href="#getting-your-alignment-objects-as-formatted-strings">6.2.2  Getting your alignment objects as formatted strings</a></li>
</ul>
</li>
<li><a class="reference internal" href="#manipulating-alignments">6.3  Manipulating Alignments</a><ul>
<li><a class="reference internal" href="#slicing-alignments">6.3.1  Slicing alignments</a></li>
<li><a class="reference internal" href="#alignments-as-arrays">6.3.2  Alignments as arrays</a></li>
</ul>
</li>
<li><a class="reference internal" href="#alignment-tools">6.4  Alignment Tools</a><ul>
<li><a class="reference internal" href="#clustalw">6.4.1  ClustalW</a></li>
<li><a class="reference internal" href="#muscle">6.4.2  MUSCLE</a></li>
<li><a class="reference internal" href="#muscle-using-stdout">6.4.3  MUSCLE using stdout</a></li>
<li><a class="reference internal" href="#muscle-using-stdin-and-stdout">6.4.4  MUSCLE using stdin and stdout</a></li>
<li><a class="reference internal" href="#emboss-needle-and-water">6.4.5  EMBOSS needle and water</a></li>
</ul>
</li>
</ul>
</li>
</ul>

  <h4>Previous topic</h4>
  <p class="topless"><a href="chr5.html"
                        title="previous chapter">Chapter 5  Sequence Input/Output</a></p>
  <h4>Next topic</h4>
  <p class="topless"><a href="chr7.html"
                        title="next chapter">Chapter 7  BLAST</a></p>
  <h3>This Page</h3>
  <ul class="this-page-menu">
    <li><a href="_sources/chr6.txt"
           rel="nofollow">Show Source</a></li>
  </ul>
<div id="searchbox" style="display: none">
  <h3>Quick search</h3>
    <form class="search" action="search.html" method="get">
      <input type="text" name="q" />
      <input type="submit" value="Go" />
      <input type="hidden" name="check_keywords" value="yes" />
      <input type="hidden" name="area" value="default" />
    </form>
    <p class="searchtip" style="font-size: 90%">
    Enter search terms or a module, class or function name.
    </p>
</div>
<script type="text/javascript">$('#searchbox').show(0);</script>
        </div>
      </div>
      <div class="clearer"></div>
    </div>
    <div class="related">
      <h3>Navigation</h3>
      <ul>
        <li class="right" style="margin-right: 10px">
          <a href="genindex.html" title="General Index"
             >index</a></li>
        <li class="right" >
          <a href="chr7.html" title="Chapter 7 BLAST"
             >next</a> |</li>
        <li class="right" >
          <a href="chr5.html" title="Chapter 5 Sequence Input/Output"
             >previous</a> |</li>
        <li><a href="index.html">Biopython_en 1.0 documentation</a> &raquo;</li> 
      </ul>
    </div>
    <div class="footer">
        &copy; Copyright 2013, Biopython.
      Created using <a href="http://sphinx-doc.org/">Sphinx</a> 1.2b1.
    </div>
  </body>
</html>